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Abstract 

We propose a semi-analytical method of calculating the timing fluctuations in mode-locked semi¬ 
conductor lasers and apply it to study the effect of delayed coherent optical feedback on pulse timing 
jitter in these lasers. The proposed method greatly reduces computation times and therefore allows 
for the investigation of the dependence of timing fluctuations over greater parameter domains. We 
show that resonant feedback leads to a reduction in the timing jitter and that a frequency-pulling 
region forms about the main resonances, within which a timing jitter reduction is observed. The 
width of these frequency-pulling regions increases linearly with short feedback delay times. We 
derive an analytic expression for the timing jitter, which predicts a monotonous decrease in the 
timing jitter for resonant feedback of increasing delay lengths, when timing jitter effects are fully 
separated from amplitude jitter effects. For long feedback cavities the decrease in timing jitter 
scales approximately as 1/r with the increase of the feedback delay time r. 
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INTRODUCTION 


Many current and future applications require ultra-high repetition frequency light pulse 
sources |RAF07| . Among these applications most also require highly regular pulse arrival 
times. Mode-locked (ML) solid state lasers can fulfill these requirements. However, such 
devices are too expensive for large scale use. Due to this limitation extensive research has 
gone into semiconductor ML lasers. The most attractive mode-locking technique, due to its 
simplicity of production and handling, is passive mode-locking, which does not require any 
external RF modulation source. However, due to the absence of an external reference clock 
passively ML lasers exhibit relatively large fluctuations in the temporal positions of pulses 
compared with a perfectly periodic pulse train [ LINlO cj. This phenomenon is referred to as 
pulse timing jitter. Recently, it was proposed to use optical feedback to significantly reduce 
the timing jitter of passively ML lasers |SOL93( ILIN10e( IOTT12a] IOTT14b] , Other meth¬ 
ods of pulse stream stabilisation which have been investigated include hybrid mode-locking 
[FIO101IARK13] and optical injection IREBlOi 1RKB 11 j . To characterize the performance 
of such devices, with respect to the timing regularity, the timing jitter is calculated. Ex¬ 
perimentally this is done using the von Linde method, which involves integrating over the 
sidebands of the power spectrum of the laser output. However, for the numerical inves¬ 
tigation of ML lasers the von Linde method can be impractical as it is computationally 
very expensive. In this paper we therefore propose a semi-analytical method of calculat¬ 
ing the pulse timing jitter for a set of delay differential equations (DDEs) proposed earlier 
to describe passive mode-locking in semiconductor lasers IVLAO 1 IVLA04aI lVLA05j . The 
method is of general nature and can be used to estimate the variance of timing fluctuations 
in a wide range of time periodic dynamical systems described by autonomous systems of 
DDEs subject to weak additive noise. 

Theoretical analysis of the influence of noise on ML pulses propagating in a laser cavity 
was first performed by H. Haus using a master equation [ HAU93a ]. Later this technique was 
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extended by taking into account the finite carrier density relaxation rate in semiconductor 
lasers [JIAOl j. The master equation has secant-shaped ML pulses as a solution, and a 
small perturbation of this state can be studied using the linearized equation of motion. The 
perturbed pulse is described by four parameters: the perturbations of the pulse amplitude, 
phase, frequency, and timing. Using the orthogonality of the solutions of the linearized 
equation to the solutions of the adjoint homogeneous linear system, coupled first order 
differential equations of motion, driven by noise, can be written out. However, due to 
multiple simplifying assumptions underlying the Haus master equation, this approach is not 
directly applicable to the analysis of semiconductor laser devices. This is why the theoretical 
estimation of timing jitter in ML semiconductor lasers has been previously performed using 
the direct numerical simulations of travelling wave [1ZHU97I . MIJL06] and delay-differential 
equation (DDE) |OTT12al IOTT14bl I.TAIJlhal ISIM14] models. As purely computational 
approaches are time-consuming, the influence of noise on the dynamics of ML pulses has 
been studied only in limited parameter regions. In a recent paper |PIM14bj a new semi- 
analytical method to estimate timing jitter in the DDE-model IVI. AO 1 IVLA04a! :VLA05 j 
of a passively ML semiconductor laser was proposed. This method was used to study the 
effect of nonlinear phenomena such as bifurcations and bistability on timing jitter, and the 
numerical results were found to be in good qualitative agreement with experimental data. In 
this paper we consider a generalisation of the semi-analytical method to study passively ML 
lasers with multiple delayed feedback. We then use this semi-analytical method to derive a 
formula for the timing jitter for resonant feedback delay lengths. 

In Section ll[]we introduce an autonomous DDE model of a laser operating in a passive 
ML regime and describe the parameters used in our calculations. In Sec. Ill, by linearizing 
the model equations near the ML periodic solution and projecting the perturbation term 
on the neutral eigenfunctions corresponding to the time and phase shift symmetries of the 
unperturbed equations, we derive a semi-analytical expression for the variance of the pulse 
timing fluctuations 1 l.\l.(i(i . 11A1.77 . Section IV is devoted to the comparison of the results 
obtained using this expression with those of direct numerical calculations of pulse timing 
jitter, and a derivation of the dependence of the timing jitter on the feedback delay time 
in the particular case of resonant feedback. Finally, in Sec. V we conclude with a brief 
discussion of our results. [ FLU07 ] 
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FIG. 1. Schematic diagram of a two section ring cavity laser subject to optical feedback from two 
external cavities (EC). The yellow region represents the gain section, the blue region corresponds 
to the saturable absorber (SA) section and the green bar indicates the spectral filtering element. 


DDE MODEL 

We use a DDE model for a passively ML ring cavity laser subject to optical feedback from 
M external cavities, based on the model introduced in [OTT12a] , a schematic diagram of 
the model is shown in Fig. [I] for the case of two feedback cavities. This model is an extension 
of the DDE model proposed in [ VLA041IVLA05] . A detailed description and derivation of 
the feedback terms for a laser with a single feedback cavity can be found in OTT 12 a . The 
final set of three coupled delay differential equations is 

£ (t) = - (7 + iu) £ (t) +7 R(t- T) e ~ i ( An +^ T £ (t 
+ 7 Em=l E *1 K m ,ie- ilC -R C t-T- Wm) e -*( An+w) (T+lTm ) £ (f — T 
G (; t ) = J g - 7 g G (; t ) - e-W (e G W - l) \£ (t) | 2 , 

Q (t) = J q - 7 q Q (t) - r s e~ Q W - l) \£ (■ t ) | 2 , 

with 

R(t) = . 

The dynamical variables are the slowly varying electric field amplitude £, the saturable gain 
G and the saturable loss Q. The saturable gain G and saturable loss Q are related to the 
carrier inversion in the gain and absorber sections, respectively. In Eq. ([2]) J g is related 


T) 

lT m ) + DC, (t) , (1) 
( 2 ) 
( 3 ) 


( 4 ) 
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to the current pumped into the gain section and J q in Eq. (J3| describes the unsaturated 
absorption. The carrier lifetimes in the gain and absorber sections are given by 1 / 7 S and 
I /7 q , respectively. The factor r s is the ratio of the saturation intensities in the gain and 
absorber sections. The M + 1 delay times in this system are the cold cavity round-trip time 
T and the external cavity round-trip times (delay times) r m of the M feedback cavites. The 
cold cavity round trip time is defined as T = v/L, where L is the length of the ring cavity. 
The bandwidth of the laser is limited by the finite width of the gain spectrum, which is taken 
into account by a Lorentzian-shaped filter function of width 7 . to describes the shift between 
the reference frequency and the central frequency of the spectral filter. The possibility of 
detuning between this latter frequency and the frequency of the nearest cavity mode is 
allowed for by the inclusion of AQ. The optical feedback is described by the sum in Eq. ([Tj). 
Here l is the number of round-trips in the external cavity, K m j is the round-trip dependent 
feedback strength of the mth feedback cavity and C m is the phase shift that accumulates 
over one round-trip in the external cavity. Below we consider feedback contributions only 
from light that has made one round-trip in the external cavities (K m , 1 = K m ). The last 
term in Eq. ([TJ) models spontaneous emission noise using a complex Gaussian white noise 
term £(f) = £1 (t) + i&it) with strength D , 

(&(*)) = 0 and (£,i(t)£,j(t')) = 5ij5(t - t'). 

Equation ([d]) describes the amplification and losses of the electric held during one round-trip 
in the laser cavity. Internal and out-coupling losses are taken into account in the attenuation 
factor k and the linewidth enhancement factors (a-factor) in the gain and absorber sections 
are denoted a g and a q , respectively. 


symbol value symbol value 


T 

25 ps 

7 

2.66 ps - 

7 g 

1 ns -1 

7g 

75 ns -1 

J 9 

0.12 ps” 1 

J, 

0.3 ps -1 

r s 

25.0 

Cm 

0 

K 

0.1 

AQ 

0 


TABLE I. Parameter values used in numerical simulations, unless stated otherwise. 
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PERTURBATION ANALYSIS 


Various methods of calculating the timing jitter are discussed in [ MIJL061 IOTT14bl 
IKEF08 , [PlMl4bl ILIN 86 ] . In this section, we consider an extension of the semi-analytical 
method of timing jitter estimation proposed in |PIM14bj for the DDE model of passively ML 
laser to the system 0-(§ with external feedback and, hence, multiple delay times. Details 
of the derivation of the semi-analytical expression for the estimation of pulse timing jitter 
are presented in the Appendix. As we do not use a specific form of equations (JT|)- (J3]) , the 
same approach can be applied to the analysis of the effect of small additive noise on stable 
periodic solutions in other physical systems described by autonomous DDEs with multiple 
delays. The advantage of the proposed method, compared with the von Linde technique 
or the so called long term jitter calculation |OTT14bj . is that it is based on the numerical 
solution of deterministic equations and therefore requires much shorter computation times. 
Furthermore, when the spontaneous emission noise is modeled by a Gaussian white noise 
term, the fluctuations of the pulse arrival times behave like a random walk |OTT14bj . making 
the timing jitter calculated from the semi-analytical method proportional to the rms timing 
jitter given by the von Linde method. This is useful for comparison with experiments. 

We consider a periodic ML solution, -0o = (Re£ 0 , Im£ 0 , G 0 , Qo) T °f the system (Jl])-([3]) 
for D = 0, with period T 0 . One should note that due to the rotational symmetry, there is 
a family of such solutions ■ -0 O = (Re(e* ¥, ^o), Ym{e ltp £o), G 0 , Qo) T , where T^, denotes the 
corresponding matrix of rotation of the So plane. The noise perturbation is assumed to 
be reasonably small, h <C 1 , and we restrict our analysis to the situation when solutions 
remain at a distance of order D from the torus of stable periodic solutions T^ ■ ij; 0 (t + 9) 
at all times (that is, the probability of a large fluctuation of the solution is assumed to be 
negligible during the typical time interval of system observation). Under this assumption, 
the noise results in a slow diffusion of the phase 9 of the solution, as well as a slow diffusion 
of the angular variable p. Furthermore, one expects that the variance of the phase 9 and 
of the variable p increases linearly with time, that is (9 — 9) 2 oc t, which expresses a simple 
diffusion process [DAF98] . We use the coefficient of proportionality in this relationship as a 
measure of the timing jitter. 

The phase of a solution can be defined in several ways |RIC02| . which, in practice, lead 
to equivalent or close results when applied for the evaluation of the phase diffusion rate. In 
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particular, the definition of the asymptotic phase is based on the fact that every solution 
of the unperturbed system 0-0 with D — 0 converges to a periodic solution P^, • -0 o (f + 9) 
in the limit t —> oo where the constant 9, called the asymptotic phase, and the angle 
ip are specific to the initial state of the solution Recall that states of system ([Tj)- 

(j3| are functions defined on the interval [—r' M , 0]. The asymptotic phase 9 and the angle ip 
remain constant along the trajectories of the unperturbed system. However, in the perturbed 
system, the asymptotic phase 9 and the angular variable ip evolve as functions of the evolving 
state t/j(t + r) (r 6 [— 

As dynamics are restricted to a small neighborhood of the limit cycle if ) 0 (and its rotations 


r v • ipo), the evolution of the phase can be deduced from the linearization (15) of system 
Eqs. <Jl])-((3]) around this cycle. Details on the analysis of the dynamics of the solutions of 
the linear system (15) as well as its effect on the evolution of the phase can be found in 
Appendix. Noise results in a slow diffusion of the variables 9 and ip along the neutral periodic 


eigenmodes of the linearized unperturbed system (16) with the variance proportional to time. 
There are two such neutral modes, 


Hi(t) = (Re£ 0 (t),lm£ 0 it),G 0 (t),Q 0 (t)) T , H 2 (t) = (-Im£ 0 (t), Re£ 0 (f), 0, 0) T , (5) 


which correspond to the time-shift and rotational symmetries of the unperturbed (D = 
0 ) nonlinear system 0-0, respectively; all the other Floquet modes are exponentially 
decaying. Two properly normalized ( 22 ) neutral modes Sijj\(t) and H\{t) of the adjoint 
linear system 0 can be used for calculating the projections of noise onto the eigendirections 
Sipi and H'i- Using the perturbation expansion with respect to the small parameter D } and 
adapting the asymptotic analysis from [ REBll J. we obtain the following equations for the 
noise-driven slow evolution of the phase 9 and the angular variable ip of solutions to Eqs. ([Tj)- 


9 = D 5i/j\(t + 9)T_ ip w(t), cp — D + 9)T- V> w(t) ( 6 ) 

with the Langevin term T_ v w(t) = (£i (t) cos ip + £ 2 (£) sin<^, —£i(t) sin ip + £ 2 (£) cos 99 , 0, 0) T 
and the T 0 -periodic coefficients and S^l- 

The coefficients of the Fokker-Planck equation for the joint probability density p(t,9,ip) 
of the stochastic process ([ 6 ]) are also periodic with respect to time. Since, for fi < 1 , the 
probability density function p(t, 9, ip) changes slowly, Eqs. (j 6 ]) and the corresponding Fokker- 
Planck equation can be averaged over the period T 0 of the functions + 9), resulting in 
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the diffusion equation with constant coefficients |KRQ91j . The diffusion coefficient 


= £ [ M,i(s )) 2 + M, 2 (s)) 2 d3 (7) 

Jo 

of the time averaged Fokker-Planck equation approximates the rate of diffusion of the phase 
9 (see Appendix). Finally, since the pulse timing jitter is usually calculated over a long time 
interval nT 0 with n 1 and the average period T 0 fa T 0 , and is normalized by the number 
of round-trips n, we make the estimate of timing jitter as the product of the diffusion rate 
by the period 

°var = dnTo = D 2 f (cb/>| a (s )) 2 + (5^[ 2 (s)) 2 ds. ( 8 ) 

Jo 

This value is approximately equal to the variance of 9(nT 0 ) divided by n 1. We note that 
for the number of roundtrips n ^ 1 that is not sufficiently large, the numerically calculated 
timing jitter is not approximated by (J 8 | since the numerically calculated value is affected by 
amplitude noise, or, in other words, stable eigendirections play role as well (see Fig. |(a)). 

For the case of resonant optical feedback, expression ([ 8 ]) for the timing jitter can be 
further simplified, to ascertain the dependence on the feedback delay length. This will be 
shown in the next section where we compare the analytic result with a numerical estimate 
of the timing jitter. 


RESULTS 

Comparison of semi-analytical and numerical methods of timing jitter calculation. 

In this section we compare the timing jitter calculated using Eq. (J8]) with that obtained 
from the variance of the pulse timing fluctuations (long-term timing jitter) through numerical 
integration of the stochastic system (Eqs. ([l|-(|3| with D ^0). The latter (numerical) 
method is described in detail in [OTT14b] . We will focus mainly on the case of one feedback 
cavity, M — 1, and compare the two approaches to the timing jitter calculation at different 
feedback delay times (r\ = r) and the feedback strengths (K { = K). 

First, we apply the semi-analytical method of the timing jitter calculation to the case 
of a passively ML semiconductor laser without feedback, i.e. K m = 0 in Eqs. ([l])-(|3|. In 
[OTT14b] it was shown that after a sufficiently large number of roundtrips n within the 
laser cavity the variance of the pulse timing fluctuations grows linearly with the round- 
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FIG. 2. (a) Comparison of the results of numerical calculation of pulse timing jitter (green solid 
line), obtained for different numbers of round-trips n, with the timing jitter value from formula ([ 8 ]) 
(red dashed line), (b) Estimation of timing jitter, calculated using formula ([ 8 ]) (red solid line) and 
the numerical method (green dots), vs noise strength D. Parameters: K = 0, r = 0, T = 25 ps, 
k = 0.3, 7 _1 = 125 fs, 'y g 1 = 500 ps, 7" 1 = 5 ps, s = 10, = 10 ps, g^ 1 = 250 ps, a g = 2, 

OLq = 1 . 


trip number. In the numerical method the timing fluctuations are therefore calculated over 
many thousands of cavity roundtrips. I11 Fig. [2] (a) the timing jitter is plotted as a function 
of the round trip number n. The initial decrease of the numerically calculated timing 
fluctuation variance (green line) with n (for small n) can be attributed to the impact of the 
eigenfunctions with Re A < 0 (see Appendix). Using DDE-BIFTOOL f ENGOl j. for 7 T 1 
(or 7r m 1), one can typically observe that many characteristic exponents A of the ML 


solution have real parts close to 0, and, therefore, the equation of motion (21) suggests that 
such exponents will have a non-negligible impact on the numerically calculated timing jitter 
even after many cavity round-trips. Since the eigenfunctions with Re A < 0 are neglected in 
the semi-analytical approach, the value of the timing jitter estimated using this approach 
does not depend on n (dashed red line in Fig. [2] (a)). In the limit of large n this value 
is in agreement with the data obtained by direct numerical integration of Eqs. ([l])-(j3]), as 
shown in Fig. [2] (a). Figure [2] (b) shows the timing jitter, obtained using both methods, in 
dependence of the noise strength D. It is seen that good quantitative agreement is obtained 
for small to moderate levels of noise. 

Next, let us consider a system with feedback from one external cavity. Figure [3] (a) 
shows a comparison of the timing jitter calculated from the two methods in dependence 
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FIG. 3. (a) Timing jitter in dependence of the noise strength, calculated using the semi-analytic 
method (red line) and the numerical method described in (OTT14b| (green dots) for r = 70Tjsi.o- 
(b) Timing jitter in dependence of the feedback delay time, calculated using the semi-analytic 
method (red dashed line) and the numerical method (green line) for D = 0.2. Parameters: a g = 0, 
a q = 0, K = 0.1. Other parameters are as in Table [TJ 


of the noise strength. For the numerical timing jitter calculation method (green dots) the 
timing fluctuations that arise over 40000 round-trips in the laser cavity are calculated, and 
the variance of these timing fluctuations is then calculated for 300 noise realisations. For 
the semi-analytical method (red line) the solutions to the adjoint linearized homogeneous 


system (17) are numerically calculated. In both cases we simulate for a sufficiently long time 
(approximately 5000 roundtrips) before starting the calculation of the timing jitter to avoid 
transient effects. We find very good agreement between the results obtained using the two 
methods. For the simulations presented in Fig. [3] (a) the feedback delay time was chosen to 
be resonant with the ML pulse repetition period (inter-spike interval time) Tjsifl of a solitary 
laser (ML laser without feedback), meaning that the condition r = qTjsi,o is fulfilled, where q 
is an integer. Resonant feedback applied in the fundamental ML regime does not significantly 
affect the dynamical behaviour of the system, hence the laser output remains periodic and 
the semi-analytic method is applicable. When the feedback delay time is tuned from one 
resonance to the next, bifurcations can occur and the dynamical behaviour can change. This 
is described in detail in [OTT14] and |OTT12a] . In Fig. [3] (b) the numerically calculated 
dependence of the timing jitter on the delay time r is compared to that estimated semi- 
analytically, spanning from the 67th to the 68th resonance (q = 67 and q = 68, respectively). 
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Feedback delay time r/T ISIfi Feedback delay time t/T IS i , 0 


FIG. 4. Timing jitter in dependence of the feedback cavity delay time and feedback strength, 
calculated numerically (a) and using the semi-analytical method (b). The timing jitter is indicated 
by the colour code and a^o is the timing jitter of the solitary laser. Regions in white indicate 
a timing jitter greater then 20fs. In subplot (b) black marks the regions where the deterministic 
system has a non-periodic solution and the semi-analytical method cannot be applied. Parameters: 
D = 0.2, a g = 2, a q = 1.5, others are as in table [TJ 


Within the frequency-pulling regions of the main resonances there is very good agreement 
between the results obtained using the two methods. The frequency pulling regions are the 
r ranges about the main resonances within which there is one pulse in the laser cavity and 
the repetition rate tunes with r |OTT12a| . In Fig. [3] (b) these regions can be identified 
by the low timing jitter about the main resonances. At the edges of the frequency-pulling 
regions there is a sharp increase in the timing jitter. This very large timing jitter coincides 
with saddle-node bifurcation points of the deterministic system (Eqs. ([]])-Q with D = 0) 
(OTT14bj . At the edge of the 67th resonance there is a large discrepancy between the semi- 
analytical and numerical methods. This is because in the stochastic system noise induced 
switching between bistable solutions, which arise due to the saddle-node bifurcations, occurs. 
Away from the bifurcation points there is good agreement between the two methods, also 
between the main resonances, because although the dynamical behaviour changes between 
the main resonances, i.e. multiple feedback induced pulses, the solutions remain periodic 
and therefore the semi-analytical method is applicable. 

For the parameters used in Fig. [3] (b) the system is well behaved and the solutions are 
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periodic, however for other parameters, particularly for larger feedback strengths and non¬ 
zero amplitude-phase coupling, this is not the case; quasi-periodic or chaotic dynamics can 
be observed. In such regions the semi-analytic approach is invalid, however the timing jitter 
calculated by numerical methods is not meaningful in these non-periodic region either. In 
Fig. [4] the timing jitter, calculated from the numerical (a) and semi-analytical (b) methods, 
is plotted in dependence of K and r for a g = 2 and a q = 1.5. The timing jitter is given by 
the colour code, where blue regions indicate a reduction in the timing jitter with respect to 
the solitary laser, red tones indicate an increase and white regions indicate a timing jitter 
greater than 20fs, indicative of a non-periodic pulse stream. In the black regions in Fig. [4] 
(b) the solutions of the DDE system are non-periodic and the semi-analytic method is not 
applied. Good agreement is observed between these two methods over most of the parameter 
range depicted. The non-periodic regions indicated in subplot (b) coincide with the very 
high timing jitter estimations obtained using the numerical method. 

A key difference between the two methods is that the semi-analytic method is based on 
the numerical simulation of deterministic equations, while the purely numerical method re¬ 
quires integration of a system of stochastic DDEs. Using the latter method one can run into 
problems that arise due to the multiplicity of stable solutions found in this system. Since 
timing jitter estimation requires averaging over many noise realisations, depending on the 
particular realisation, due to transient effects, the system can land on different solutions. As 
different ML solutions can have slightly different inter-spike interval times, the fully numer¬ 
ical estimation of the timing jitter can lead to erroneously large values in such case [ SIM14 J. 
This makes it difficult to perform timing jitter calculations over a large parameter domain, 
as it is not easy to distinguish between the above mentioned effect and a destabilisation 
of the pulse stream due to the feedback conditions. Note that this is a different effect to 
switching between solutions within one time series. Such difficulties are eliminated when 
using the semi-analytic method, as in this case the estimation of the variance is based on the 
integration of deterministic equations. Therefore, there are two main advantages to using the 
semi-analytic method to calculate the timing jitter, compared with brute force methods in¬ 
volving numerical integration of stochastic differential equations. Firstly, the aforementioned 
difficulties can be avoided, and secondly, the computation times can be greatly reduced (by 
over a factor of 100) as averaging over many noise realisation is not needed. This means that 
it can become feasible to calculate the timing jitter for longer feedback delay times, which 
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is of interest due to the improved timing jitter reduction predicted for increased delay times 
[OTT14] and for better comparison with experiments, where typically very long feedback 
cavities are used | ARSl3l iLlNlOe J. 


Delay length dependence of timing jitter 
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FIG. 5. (a) and (b) Timing jitter an in dependence of the feedback cavity delay time. The colour 
code indicates the timing jitter according to the colour bar given in subplot (c). The black dashed 
line indicates the timing jitter of the solitary laser, (c) Timing jitter an in dependence of the 
feedback cavity delay time, where r = tq + n for any given point. The horizontal axis spans one 
Tisi,t= o an d is centered on an exact main resonance. The vertical axis indicates the number of 
the main resonance. The timing jitter is indicated by the colour code and an o is the timing jitter 
of the solitary laser. Parameters: K = 0.1, D = 0.2, a g = 0, a q = 0, others as in table [TJ 


We now use the semi-analytic method to investigate how the timing jitter decreases with 
increased resonant feedback delay times and how the width of the frequency-pulling regions 
is affected by this increase. In Fig. [5] the timing jitter is plotted as a function of r in subplots 
(a) and (b) for a short and a long r range, respectively. The black dashed line indicates 
the timing jitter of the solitary laser. The delay times are plotted in units of Tjsi, t = o , the 
inter-spike interval time for zero delay feedback (instantaneous feedback, r = 0 and K ^ 0), 
meaning that the resonant feedback occurs at the integer delay values. ( Tjsi, t =o and Tjsi, o 
only differ slightly. Here we choose Tjsi, t =o as our reference because the period is the same 
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for all r = qTjsi, r = o, g £ N, and we will use this property in subsequent calculations.) In 
both (a) and (b) a timing jitter reduction is observed for resonant feedback. For the longer 
delay times depicted in subplot (b) the timing jitter reduction is greater and the frequency¬ 
pulling region about the main resonances is wider. Changes in the frequency-pulling regions 
are not discernible over small r ranges. To show the change in dependence of r more clearly 
a map of the timing jitter is shown in a r — r plot in subplot (c). In this plot both axes 
are related to the delay time, the Ti axis shows changes over one T/5/ iT=0 -interval , whereas 
the t 0 axis shows changes from one resonance to the next. For each point on this map the 
feedback delay time is given by r = To + T\. The T\ axis is centered on the exact main 
resonances r = qTisi, r = o and the To axis gives the number q of the main resonances. The 
timing jitter is given by the colour code. Regions in blue and green indicate a reduction in 
the timing jitter with respect to the solitary laser (K = 0) and regions in red indicate an 
increase in the timing jitter. In the green regions the timing jitter is reduced by a factor 
of 10 or greater. For all q values a reduction in the timing jitter is achieved at the exact 
main resonances and for increasing q the decrease in the timing jitter can clearly be seen. 
It is seen from Fig. |5]^c) that for short delays the width of the frequency-pulling regions, 
with reduced timing jitter, increases approximately linearly with the number q. The edges 
of the frequency pulling region are marked by the dashed black lines. At about q = 50 
the frequency-pulling region is intersected by the solutions that correspond to higher order 
resonances ( pr = qT IS p r = o, where p = 2,3,4,...). This is due to a bistability between the 
main and higher order resonant solutions |OTT12aj . For the results presented in subplot 
(c) of Fig. [5j the same initial conditions were used in the numerical simulations for all delay 
values. By performing a sweep in r (using the previous r solution as the initial conditions 
for the next r) one can stay on the main resonant solution in the bistable regions. 


In order to quantify the decrease in the timing jitter with increasing number q , we have 
plotted the timing jitter at the main resonances in Fig. [6j The red line shows the results of the 
semi-analytic method for the exact main resonances r = qTjsi , r =o ( r values corresponding 
to the white dashed line in Fig. [5] (b)) and the blue line shows the results of the semi- 
analytic method for the minimum timing jitter in each main resonance frequency-pulling 
region (r values corresponding to the white dot-dashed line in Fig. [5] (b)). The expression 
for the timing jitter at the main resonances, r = qTjsi,r= o, can be derived analytically using 
Eq. Q and the bilinear form (18). At the exact main resonances the solutions to Eqs. ([!])- 
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FIG. 6. Timing jitter ait at the exact main resonances (red solid line) and the minimum timing 
jitter in each resonance region (blue solid line) as a function of the number q of the main resonance, 
calculated using the semi-analytic method. The dashed line shows the timing jitter at the exact 


main resonances given by the analytic expression Eq. (13). The dot-dashed line shows the fit of 


Eq. (14) to the minimum timing jitter in each resonance regions. Parameters: K = 0.1, D = 0.2, 


a g = 0, a q = 0, others as in table [TJ 


(J3| are identical for all q, and the periodicity is the same as that of the laser with zero 


delay (instantaneous) feedback T 0 (Tisi, t = o)- Therefore, for r = qTjsi, T =o, E<1- (18) can be 
expressed as 

r 0 

\Sijj\ S'ljj] ( t ) = Sijj^{t)Sxjj{t) + / Sift (t + r + T) B 0 (t + r ) S'lp (t + r) dr 

r° 

+K / §^ (t + r + T) Bi (t + r ) d'lp (t + r) dr 


' —T 


+K f 5^ (t + r + T ) B 1 ( t + r) 5^ ( t + r) dr. 

J -T-qTiSi, t=o 


( 9 ) 


The last term on the right-hand side can be further simplified due to the time shift invariance 
and periodicity of the integrand, giving 

c0 


= \5ijj\ Sijj] T +Iiq {Sijj^ (t + r + T)) B x (t + r) S'ljj (t + r) dr , (10) 

•J —TiSI,t=0 


where the first three terms on the rhs of Eq. ([9]) are now expressed as [Sijj\ Sijj] T °, which 
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is the bilinear form for r = 0 (q — 0). Equation (j8l) can thus be expressed as 


&vn.r 


/•To 

D 2 / 

Jo 

will (t) 

2 

+ 

(/) ) 


" ° + KqT'(K) J 

V 


T ° + KqT’ (K) j 


dt, 

( 11 ) 


where \WM] = 5 ^ ] and 


r (. k) = f 

J-: 


(t + r + T) B 1 (t + r) 5ip * (t + r) dr, 

' —TlSI,T =0 

which is a function of K but not of r. Finally, Eq. 0 can be simplified to 


rT 0 


@nn.r 


1 + KqT (. K) 


D 2 


foot,! 0 ®) +(^oI,2°W) dt, (12) 


where 8ij )° = 

ity condition for r = 0 and T (K) = -— ^ A ^ T=0 . 


^i 0 , ^[,2°, S4 T t/, ^V ; oM°j i s l^ ie solution fulfilling the biorthogonal- 

r . The timing jitter for resonant feedback, 

Y>w{ ,°n> ij 

r = qTjsi )T = o, is therefore given by 


cr 


T = g7/SI,r=0 


O', 


T—0 


( A ') 


(13) 


where cr^ 0 


1 + KqT {K) ’ 

(Jl) is the timing jitter for r = 0. The curve obtained using this analytic 
expression is shown by dashed black line in Fig. [6j A formula for the minimum jitter can 
not be derived in the same way as the inter-spike interval time changes with q. However, 
fitting the minimum jitter curve for various feedback strengths we find that the relation 


Tt 


cr, 


T—0 


(K) 


1 + Kq’ (14) 

holds well for low feedback strengths. The fit is plotted in the black dot-dashed line in Fig. [6j 


In the derivation of Eq. (13) contributions to the timing jitter from eigenfunctions with 


negative eigenvalues, A < 0, are neglected. However, for increased feedback delay lengths, 
the number of weakly stable Floquet multipliers close to one increases. This leads to long 
transients in numerical simulations of the deterministic system (Eqs. Q-Q D — 0). These 
transient effects are accompanied by fluctuations in the pulse heights, which have the pe¬ 
riodicity of the feedback delay time. Including noise in the system excites these transient 
amplitude fluctuations, which results in an increased timing jitter, as, via the interaction 
with the gain and absorber media, changes in the pulse height also lead to slight changes in 
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the pulse positions. Equation (13) is therefore only valid in the limit in which such effects 


can be neglected. For the parameter values used in our simulations Eq. (13) holds for up to 
q ~ 300. These noise induced transient effects were observed experimentally as side peaks 
in the phase noise spectra [ HA.T12 . 1ARS131 DRZ13a j. 


CONCLUSIONS 


We have investigated the influence of optical feedback on the timing jitter of a passively 


ML semiconductor laser. For resonant feedback we have derived an expression, Eq. (13), for 


the analytical dependence of the timing jitter on the feedback delay length, showing that the 
timing jitter drops off as approximately 1/r for r T, as long as amplitude jitter effects 
can be neglected. About the main resonant feedback delay lengths, frequency-pulling re¬ 
gions form, in which the timing jitter is reduced with respect to the solitary laser. For small 
feedback strengths K the widths of these frequency-pulling regions increase linearly with 
the number q of the main resonance. These results were obtained using a semi-analytical 
method, presented in this paper, of calculating timing fluctuations in a DDE system describ¬ 
ing the dynamics of a passively ML semiconductor laser subject to optical feedback from an 
arbitrary number of feedback cavities. The semi-analytical method shows good agreement 
with methods based on direct numerical integration of the stochastic model and has the 
advantage of greatly reduced computation times. 
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Appendix: Derivation of the expression for the rate of the phase diffusion 

Here we derive formula (Jr]) for the phase diffusion rate. Recall that V’o(^) is a To-periodic 
ML solution of system @-§f Substituting the expression -0( t ) = 'ipo (t) + into this 
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system, we obtain the linearized equations 


M 


—5^ (t) = A (t) Sip (t) + > v B m (t - r' m ) Sip (t - r' m ) + Dw(t), 


m=0 


(15) 


where A and B m are T 0 -periodic Jacobi matrices of the linearization; Tq = T , T' m = T + r m 
for m ^ 1; and, Dw(t ) = D(£i(t), £ 2 (t), 0, 0) T is the small noise term. Explicit expressions 
for the matrices Apt) and B(t ) can be found in [ REBll J. When there is no noise (D — 0), 
the homogeneous system 

, M 

- Bm - r m)<W - T 'm) = 0 (!6) 


and its adjoint system, for a row vector Sip^pt) = ( Sip \, Sipl, Sip' 3 , Sip\), 


d 

dt 


M 


5^(t) + 5^{t)A(t) + y + 4)S m (() = 0, 


(17) 


m=0 


have characteristic solutions (eigenmodes) of the form Sippt) = Sip\{t) = e xt p\(t ) and 
Sip\t ) = Sip\(t) = e _At pA(^)i respectively, where functions Pa(^) and p^(f) are T 0 -periodic 


and the complex value A is a Floquet exponent of (16). The bilinear form 1 1 AI.66 . HAL7 7] 


M 0 


[Sip\ Sip] (t) = Sip\t)Sip(t) + ^ + r + T' m )B m (t + r)Sip(t + r)dr (18) 


171= 1 T m 


is instrumental in quantifying the effect of noise along different eigenclirections Sip\(t) for the 


perturbed system (15), because for every solution Sippt) of (15) and every solution Sip\t) of 


(pT|) the following relation holds at all times: 


d[Sip\ Sip](t) 
dt 


= DSip^ (t)w(t). 


(19) 
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Indeed, 


(< 


^[6^,5ip](t) = I S^( s + t + r L)B m (s + t)5'ilj(s + t)ds 


dS ^dt^ S ^ + H\t) dS ^- + J t Y1 \ 5 ^( s + T , m ) B m (s)8'i/;(s)ds 

m Jt ~ T m 

= ~ + Y1 + T 'm) B ™{t) j H{t) 

+H\t) | A(t)5ip(t) + ^B m (t - r'JS^it - T' m )+w{t) 

+ ^iH\t + ^ m ) B mit)Hit) - 5^(t)B m (t - r'JS^it - t'J) 

m 

= D Si/j^ {t)w(t). 


In particular, for every pair of solutions of the homogeneous systems (16) and (17) (D = 0), 


the form [8^\ H\{t) is independent of time. Eq. ( |19j ) also ensures the biorthogonality 
property 

[H\ ,<%](*) =0 ( 20 ) 


for any pair of eigenfunctions of problems (16) and (17) with A 7 - fi. Furthermore, Eq. (19) 


implies that for any solution Sipit) of the inhomogeneous problem (15), the projection y\(t) = 
e xt [8ip\, H]{t) satisfies the equation 

dy\(t) 


dt 


= + Dp\{t)w{t) 


( 21 ) 


with the Langevin term w(t). For Re A < 0, this equation defines an Ornstein-Uhlenbeck 
type process with a uniformly bounded variance of order D 2 . On the other hand, for A = 0, 
we obtain a process similar to the Brownian motion with the variance that grows linearly with 


time as D 2 t. Hence, noise mostly affects the projections of a solution of (15) onto the neutral 


eigenmodes (|5]) that have A = 0. The two corresponding adjoint neutral eigenfunctions (that 
is, T 0 -periodic solutions of the adjoint system ©) can be normalized in such a way as to 
satisfy the relations 

Hi,Hi {t) = 5ipl,Hi it) = 1, H 11 H 2 it) = Hl,Hi (f) = 0. (22) 

For stable mode-locked solutions ipo it) all the non-zero Floquet exponents of the linearized 
system have negative real parts. 
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Using the linearization, we can approximate the asymptotic phase of a solution to the 
nonlinear system (Jl])-(J3]) by the formulas 


- $) --00 (t + 0)= 5^>' 2 ,r_^(t-d) — ip 0 (t + $) = o, 


(23) 


These equations define the “time” phase 9 and the “angular” phase ip implicitly for any given 


state ip(t + r ) (r e 0]) of the nonlinear system. Geometrically, (23) is a codimension 2 

linear subspace which is tangent to the surface of constant asymptotic phases 9, p at the point 
where this surface intersects the torus of shifted periodic solutions ■ ip 0 (t + 9) in the state 
space of the system. As we consider solutions that remain within a small distance of order D 


from this torus, the error between the asymptotic phase and its approximation (23) is of next 
order D 2 . Also, note that Eqs. (23) themselves can be used as an alternative definition of the 


phase, because these equations define a foliation of a small tubular neighborhood surrounding 
the torus of periodic solutions by non-intersecting surfaces 9 = const, p = const. 


In order to derive the equation for the evolution of the phase, we differentiate Eqs. (23) 


with respect to t,9,p. Using symmetry, one obtains from Eq. (19) the relationship 
d 


dt 


Si/)}, Y_ v ip(t - 9) - ip o (t + 9) = D Sip}(t + 9)Y_ ip w(t) (24) 

for i — 1,2. When differentiating the bilinear form Sip}, Y^ipft — 9) — ipo (t + 9) with 
respect to 9 and ip, we omit the terms that are proportional to ip — Y v ip 0 (t — 9), because 
these therms have the order D in the small vicinity of the cycle that we consider. In this 
approximation, we obtain 


d_ r 
89 

d_ r 

dtp 


HlJ-v'Kt-o) — -0o 


Srpl,r_ v i’(t -e) 


0 t + 9) — — 

Sip}, Sipx 

{t - 9) - — 

Slp\,Slp 2 


(25) 


(26) 


Combining relationships (22)—(26), we arrive at the coupled system of stochastic equations 


i§ that describe the slow evolution of the variables 9 and <p. 

Finally, using the Feynman-Kac formula, we obtain the Fokker-Planck equation for the 
joint probability density p(t, 9, ip) of the stochastic process (§: 


dp fid 2 

dt ~ V 2 m 2 ^ 11 ^ + d9dp 


d 2 1 d 2 

(d 12 p) + --r^{.d 2 2 P) 
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This equation has variable diffusion coefficients 


* 1 = £> 2 (M,i) 2 + M.2) 2 )(i + «). 
d 2 . 2 = D 2 ((H{if + (^{ 2 ) 2 )(t + 6), 

<*i2 = d 3 (sii\ + s^i\ 2 Sipl 2 )(t + 0 ), 


where Sfr, are the coordinates of the 4-dimensional vector-functions Since, for D <C 1, 


the probability density changes slowly, Eq. (27) can be averaged over the period T 0 of the 
functions dij(t + 9), resulting in the diffusion equation with constant coefficients di 3 (see, for 
example, (KROQl j). The averaged coefficient d\ i that approximates the rate of diffusion of 
the phase 6 is defined by formula ((Tj) . 
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